#pragma once
#include<fstream>
#include<iostream>
#include<string>
#include<bits/stdc++.h>
#include<cmath>
#include <seal/seal.h>
#include <algorithm>
#include <sys/resource.h>
#include <sys/wait.h>
#include <unistd.h>

// #define PI 3.14159265358979323846;
// #define EARCH_R 6731000;//地球半径

using namespace std;
using namespace seal;

// const double PI=3.14159265358979323846;
const double PI=M_PI;
const double EARCH_R=6731e3;//米

const double eps = 1e-5;//用于判断两个浮点数是否相等
// const double eps = 1;
bool isEqual(double a, double b) {
    return fabs(a - b) < eps;
}

//定义工作者和任务
enum WT_type{
    worker,//0
    task,
};
//工作者和任务花费
union WT_cost
{
    double pay;//支付
    double reward;//获得的报酬
    double rate;
};
vector<pair<double,double>> Xposition;//存加密后位置，用于计算相关矩阵
vector<vector<double>> SMatrix(2,vector<double>(2,0));//协方差矩阵

//节点定义
struct Node{
    WT_type type;
    pair<double,double> loc;
    double radius;
    int capacity;
    int beginTime,endTime;
    WT_cost cost;
    void show(){
        if(type==worker){
            cout<<beginTime<<"\t"<<type<<"\t"<<loc.first<<"\t"<<loc.second<<"\t"<<radius<<"\t"<<capacity<<"\t"<<endTime<<"\t"<<cost.rate<<endl;
        }else{
            cout<<beginTime<<"\t"<<type<<"\t"<<loc.first<<"\t"<<loc.second<<"\t"<<endTime<<"\t"<<cost.pay<<endl;
        }
    }
};


/**
 * @brief 计算两个坐标之间的距离平方
 * 
 * @return double 
 */
inline double calDistance(pair<double,double>loc1,pair<double,double>loc2){
    return (loc1.first-loc2.first)*(loc1.first-loc2.first)+(loc1.second-loc2.second)*(loc1.second-loc2.second);
}
/**
 * @brief 计算两个坐标之间的距离
 * 
 * @return double 
 */
inline double calDistance2(pair<double,double>loc1,pair<double,double>loc2){
    return sqrt((loc1.first-loc2.first)*(loc1.first-loc2.first)+(loc1.second-loc2.second)*(loc1.second-loc2.second));
}
inline double calDistance3(pair<double,double>loc1,pair<double,double>loc2){
    // double avg_x = 0,avg_y = 0,covxx = 0,covxy = 0,covyx = 0, covyy = 0;
    // avg_x = (loc1.first-loc2.first)/2;
    // avg_y = (loc1.second-loc2.second)/2;
    // covxx = pow(loc1.first-avg_x,2)+pow(loc2.first-avg_x,2);
    // covxx /= 2; 
    // covyy = pow(loc1.second-avg_y,2)+pow(loc2.second-avg_y,2);
    // covyy /= 2; 
    // covxy = (loc1.first-avg_x)*(loc1.second-avg_y) + (loc2.first-avg_x)*(loc2.second-avg_y);
    // covyx = covxy;
    double dist,a1,a2;
    a1 = loc1.first-loc2.first;
    a2 = loc1.second-loc2.second;
    dist = pow(a1,2)*SMatrix[0][0]+a1*a2*SMatrix[0][1]+a1*a2*SMatrix[1][0]+pow(a2,2)*SMatrix[1][1];
    // dist = pow(a1,2)*covxx+a1*a2*covxy+a1*a2*covyx+pow(a2,2)*covyy;
    // dist = sqrt(dist);
    return dist;
}

/**
 * @brief 检查是否满足在半径内
 * 
 * @param worker 
 * @param task 
 * @return true 
 * @return false 
 */
inline bool satisfyDistance(Node& worker,Node& task){
    if (calDistance(worker.loc,task.loc) > worker.radius*worker.radius)
    {
        return false;
    }
    return true;
}
inline bool satisfyDistance2(Node& worker,Node& task){
    double dist,a1,a2;
    a1 = worker.radius;
    a2 = worker.radius;
    dist = pow(a1,2)*SMatrix[0][0]+a1*a2*SMatrix[0][1]+a1*a2*SMatrix[1][0]+pow(a2,2)*SMatrix[1][1];
    // dist = sqrt(dist);
    if (calDistance3(worker.loc,task.loc) > dist)
    {
        return false;
    }
    return true;
}
/**
 * @brief 判断工作者和任务之间是否符合时间要求
 * 
 * @param worker 
 * @param task 
 * @return true 
 * @return false 
 */
inline bool satisfyTime(Node& worker,Node& task){
    if (worker.beginTime < task.endTime && task.beginTime < worker.endTime)
    {
        return true;
    }
    return false;
}
inline bool satisfyDistanceTime(Node& worker,Node& task){
    if(satisfyDistance(worker,task)&&satisfyTime(worker,task)){
        return true;
    }
    return false;
}
inline bool satisfyDistanceTime2(Node& worker,Node& task){
    if(satisfyDistance2(worker,task)&&satisfyTime(worker,task)){
        return true;
    }
    return false;
}


// 高斯-约旦消元求逆矩阵
int inverse(vector<vector<double>>& A)
{
    int n = A.size();
    double matrixEps = 1e-6,tmp;
    vector<vector<double>> B(n, vector<double>(n, 0));
    for (int i = 0; i < n; ++i) B[i][i] = 1;
    
    for (int i = 0; i < n; ++i)
    {
        //寻找第 i 列不为零的元素
		int ii;
		for ( ii= i;ii < n; ii++)
		{
			if (fabs(A[ii][i]) > matrixEps) //满足这个条件时，认为这个元素不为0
				break;
		}
        if (ii<n){
            if(ii!=i){ // 如果主元为0，则需要交换行
                for(int j=0;j<n;j++){
                    tmp = A[i][j]; A[i][j] = A[ii][j]; A[ii][j]=tmp;
                    tmp = B[i][j]; B[i][j] = B[ii][j]; B[ii][j]=tmp;
                }
            }
        }else{
            return -1;// 矩阵不可逆
        }
        
        // 将主元归一化
        tmp = A[i][i];
        for (int j = i; j < n; ++j){
            A[i][j] /= tmp;
            B[i][j] /= tmp;
        }
        for(int j=0;j<i;++j)B[i][j] /= tmp;
        
        // 消去第i列的其他元素
        for (int j = 0; j < n; ++j)
        {
            if (j == i) continue;
            tmp = A[j][i];
            for (int jj = i; jj < n; ++jj){
                A[j][jj] -= tmp * A[i][jj];
                B[j][jj] -= tmp * B[i][jj];
            }
            for(int jj=0;jj<i;++jj) B[j][jj] -= tmp * B[i][jj];
        }
    }
    A = B;
    return 1;
}

//计算马氏距离协方差矩阵
void calSMatrix(){
    int avg_n = Xposition.size();
    double avg_x = 0,avg_y = 0,covxx = 0,covxy = 0,covyx = 0, covyy = 0;
    for(int i=0;i<avg_n;i++){
        avg_x += Xposition[i].first;
        avg_y += Xposition[i].second;
    }
    avg_x /= avg_n;
    avg_y /= avg_n;
    for(int i=0;i<4;i++){
        covxx += pow(Xposition[i].first-avg_x,2);
        covyy += pow(Xposition[i].second-avg_y,2);
        covxy += (Xposition[i].first-avg_x)*(Xposition[i].second-avg_y);
    }
    covxx /= (avg_n-1);
    covyy /= (avg_n-1);
    covxy /= (avg_n-1);
    covyx = covxy;
    SMatrix[0][0] = covxx;
    SMatrix[0][1] = covxy;
    SMatrix[1][0] = covyx;
    SMatrix[1][1] = covyy;
    int re = inverse(SMatrix);
    // cout<< SMatrix[0][0] <<" "<< SMatrix[0][1]<<endl;
    // cout<< SMatrix[1][0] <<" "<< SMatrix[1][1]<<endl;
    if (re==-1) cout<<"not inverse"<<endl;
}
// 马氏距离
double MahalanobisDistance(vector<vector<double>> &X){
    vector<vector<double>> S(2,vector<double>(2,0));
    int avg_n = X.size();
    // cout<<"avgn:"<<avg_n<<endl;
    // X[0].push_back(3);
    // X[0].push_back(4);
    // X[1].push_back(5);
    // X[1].push_back(6);
    // X[2].push_back(2);
    // X[2].push_back(2);
    // X[3].push_back(8);
    // X[3].push_back(4);

    double avg_x = 0,avg_y = 0,covxx = 0,covxy = 0,covyx = 0, covyy = 0;
    for(int i=0;i<avg_n;i++){
        avg_x += X[i][0];
        avg_y += X[i][1];
    }
    avg_x /= avg_n;
    avg_y /= avg_n;
    // cout<<avg_x<<endl; 
    // cout<<avg_y<<endl; 
    for(int i=0;i<4;i++){
        covxx += pow(X[i][0]-avg_x,2);
        covyy += pow(X[i][1]-avg_y,2);
        covxy += (X[i][0]-avg_x)*(X[i][1]-avg_y);
    }
    covxx /= (avg_n-1);
    covyy /= (avg_n-1);
    covxy /= (avg_n-1);
    covyx = covxy;
    S[0][0] = covxx;
    S[0][1] = covxy;
    S[1][0] = covyx;
    S[1][1] = covyy;
    // cout<< S[0][0] <<" "<< S[0][1]<<endl;
    // cout<< S[1][0] <<" "<< S[1][1]<<endl;
    int re = inverse(S);
    // if (re==-1) cout<<"not inverse"<<endl;
    if (re==-1) return -1;
    // cout<< S[0][0] <<" "<< S[0][1]<<endl;
    // cout<< S[1][0] <<" "<< S[1][1]<<endl;
    
    double dist=0,dist1,dist2,a1,a2,b1,b2;
    a1 = X[0][0]-X[1][0];
    a2 = X[0][1]-X[1][1];
    b1 = X[2][0]-X[3][0];
    b2 = X[2][1]-X[3][1];
    // cout<<"a1:"<<a1<<" a2:"<<a2<<endl;
    dist1 = pow(a1,2)*S[0][0]+a1*a2*S[0][1]+a1*a2*S[1][0]+pow(a2,2)*S[1][1];
    dist1 = sqrt(dist1);
    dist2 = pow(b1,2)*S[0][0]+b1*b2*S[0][1]+b1*b2*S[1][0]+pow(b2,2)*S[1][1];
    dist2 = sqrt(dist2);
    dist = dist1+dist2;
    // cout<<"d1:"<<dist1<<" d2:"<<dist2<<" d:"<<dist<<endl;
    return dist;
}
double MahalanobisDistance2(pair<double,double>& worker1,pair<double,double>& worker2,pair<double,double>& task1,pair<double,double>& task2){
    double dist=-1,dist1,dist2,a1,a2,b1,b2;
    a1 = worker1.first-worker2.first;
    a2 = worker1.second-worker2.second;
    b1 = task1.first-task2.first;
    b2 = task1.second-task2.second;
    // cout<<"a1:"<<a1<<" a2:"<<a2<<endl;
    dist1 = pow(a1,2)*SMatrix[0][0]+a1*a2*SMatrix[0][1]+a1*a2*SMatrix[1][0]+pow(a2,2)*SMatrix[1][1];
    dist1 = sqrt(dist1);
    dist2 = pow(b1,2)*SMatrix[0][0]+b1*b2*SMatrix[0][1]+b1*b2*SMatrix[1][0]+pow(b2,2)*SMatrix[1][1];
    dist2 = sqrt(dist2);
    dist = dist1+dist2;
    // cout<<"d1:"<<dist1<<" d2:"<<dist2<<" d:"<<dist<<endl;
    return dist;
}


//原始的工作者和任务列表
vector<Node> workerLists,taskLists;
//扰动之后的工作者和任务列表
vector<Node> perturbWorkers,perturbTasks;

//全局最大流匹配成功的
vector<int> mWorkers;
vector<int> mTasks;
//扰动位置情况下的匹配对
vector<pair<int,int>> M0;

//扰动的水平l和半径r
int l=400,r=1000,k=2,lamda=2;

//同态加密 a服务器
void secureDistanceCalA(int w,int t, double& scale, CKKSEncoder& ckks_encoder,Encryptor& encryptor,Evaluator& evaluator,Ciphertext& cipherxResult){
    Plaintext plainx1,plainx2,plainy1,plainy2;
    double x1,x2,y1,y2;
    x1 = workerLists[w].loc.first;
    y1 = workerLists[w].loc.second;
    x2 = taskLists[t].loc.first;
    y2 = taskLists[t].loc.second;
    // cout<<"real:"<<(x1-x2)*(x1-x2)+(y1-y2)*(y1-y2)<<endl;
    ckks_encoder.encode(x1,scale,plainx1);
    ckks_encoder.encode(-x2,scale,plainx2);
    ckks_encoder.encode(y1,scale,plainy1);
    ckks_encoder.encode(-y2,scale,plainy2);
    
    Ciphertext cipherx1,cipherx2,ciphery1,ciphery2;
    encryptor.encrypt(plainx1,cipherx1);
    encryptor.encrypt(plainx2,cipherx2);
    encryptor.encrypt(plainy1,ciphery1);
    encryptor.encrypt(plainy2,ciphery2);

    evaluator.add_inplace(cipherx1,cipherx2);
    evaluator.add_inplace(ciphery1,ciphery2);
    evaluator.square_inplace(cipherx1);
    evaluator.square_inplace(ciphery1);
    evaluator.add(cipherx1,ciphery1,cipherxResult);
}

//同态加密 b服务器
bool secureDistanceCalB(int w,int t){
    EncryptionParameters parms(scheme_type::ckks);
    size_t poly_modulus_degree = 8192;
    parms.set_poly_modulus_degree(poly_modulus_degree);
    parms.set_coeff_modulus(CoeffModulus::Create(poly_modulus_degree, { 60, 40, 40, 60 }));
    SEALContext context(parms);
    KeyGenerator keygen(context);
    auto secret_key = keygen.secret_key();
    PublicKey public_key;//生成全局的公钥
    keygen.create_public_key(public_key);

    Encryptor encryptor(context, public_key);
    Decryptor decryptor(context, secret_key);
    Evaluator evaluator(context);
    CKKSEncoder ckks_encoder(context);
    double scale = pow(2.0,40);
    
    Ciphertext cipherx1;

    Plaintext disResult;
    secureDistanceCalA(w,t,scale,ckks_encoder,encryptor,evaluator,cipherx1);
    decryptor.decrypt(cipherx1,disResult);
    vector<double> dec_values;
    ckks_encoder.decode(disResult, dec_values);
    double d = dec_values[0];
    d = sqrt(d);
    // cout<<"HE:"<<d<<endl;
    if(perturbWorkers[w].radius>=d&&satisfyTime(perturbWorkers[w],perturbTasks[t])){
        return true;
    }
    return false;
}



void usage(){
    cout<<"Usage: ./main 40 1000 filePath"<<endl;
    cout<<"First parameter is Geo-I privacy level l"<<endl;
    cout<<"Second parameter is privacy within radius r"<<endl;
    cout<<"Third parameter is data path"<<endl;
}

void usage2(){
    cout<<"Usage: ./main 40 1000 2 2 filePath"<<endl;
    cout<<"First parameter is Geo-I privacy level l"<<endl;
    cout<<"Second parameter is privacy within radius r"<<endl;
    cout<<"Third parameter is group numbers k"<<endl;
    cout<<"fourth parameter is iterates lamda"<<endl;
    cout<<"Fifth parameter is data path"<<endl;
}

//测运行时间和内存
struct program_t {
	struct timespec real;
	struct timespec proc;
};
int usedMemory;


void save_time(program_t& prog) {
	clock_gettime(CLOCK_REALTIME, &prog.real);
	clock_gettime(CLOCK_PROCESS_CPUTIME_ID, &prog.proc);
}

double calc_time(const program_t& st, const program_t& ed) {
	double ret;

	#ifdef USING_PROC_TIME
	double proc_st = st.proc.tv_sec + (float)st.proc.tv_nsec / 1e9;
	double proc_ed = ed.proc.tv_sec + (float)ed.proc.tv_nsec / 1e9;
	ret = proc_ed - proc_st;
	#else
	double real_st = st.real.tv_sec + (float)st.real.tv_nsec / 1e9;
	double real_ed = ed.real.tv_sec + (float)ed.real.tv_nsec / 1e9;
	ret = real_ed - real_st;
	#endif /* USING_PROC_TIME */

	return ret;
}

int get_proc_status(int pid, const char * mark) {
	const int BUFFER_SIZE = 5120;
	FILE * pf;
	char fn[BUFFER_SIZE], buf[BUFFER_SIZE];
	int ret = 0;

	sprintf(fn, "/proc/%d/status", pid);
	pf = fopen(fn, "re");
	int m = strlen(mark);
	while (pf && fgets(buf, BUFFER_SIZE - 1, pf)) {

		buf[strlen(buf) - 1] = 0;
		if (strncmp(buf, mark, m) == 0) {
			sscanf(buf + m + 1, "%d", &ret);
		}
	}

	if (pf)
		fclose(pf);

	return ret;
}

void watchSolutionOnce(pid_t pid, int& usedMemory) {
	int tempMemory;
	int status;
	struct rusage ruse;
    
	tempMemory = get_proc_status(pid, "VmRSS:");
	if (tempMemory > usedMemory)
		usedMemory = tempMemory;
	
	wait4(pid, &status, 0, &ruse);
	tempMemory = get_proc_status(pid, "VmPeak:");
	if (tempMemory > usedMemory)
		usedMemory = tempMemory;
}